--- title: Peak slice maps (working on it) keywords: fastai sidebar: home_sidebar summary: "Slicing the cube " description: "Slicing the cube " nb_path: "notebooks/60_peakmaps.ipynb" ---
import skimage.exposure as ske
from maxrf4u import DataStack, HotmaxAtlas, get_hotslices, get_peakmaps, multi_plot, plot_peakslices, plot_ptrn
ds = DataStack('RP-T-1898-A-3689.datastack')
x_keVs = ds.read('maxrf_energies')
y_max = ds.read('maxrf_maxspectrum')
cube = ds.read('maxrf_cube', compute=False) # don't load into memory yet (too big)
imvis_reg_highres = ds.read('imvis_reg_highres')
imvis_reg = ds.read('imvis_reg')
imvis_extent = ds.read('imvis_extent')
hma = HotmaxAtlas('RP-T-1898-A-3689.datastack')
n = 5
x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n]
hot_pixel =hma.hotmax_pixels[n]
slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
ax, labels = hma.plot_spectrum(n)
plot_peakslices(x_keVs, slices, ax=ax);
ymin, ymax = ax.get_ylim()
ax.set_ylim(-5, ymax)
plot_ptrn('Ca', -2, ax);
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack', verbose=True)
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')
peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]
fig, axs = multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel)
Let's check if peak #5[6] is an escape peak for Ca_Ka by calculating the energy shift....
peak_idxs = hma.peak_idxs_list[5]
x_keVs[peak_idxs[0]] - x_keVs[peak_idxs[6]]
Indeed!!
title = axs.flatten()[6].get_title() + '(escape)'
axs.flatten()[6].set_title(title);
There is much more to be learned from this overview, but before going into this let's look at other elements.
n = 9 # Mn
x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n]
hot_pixel =hma.hotmax_pixels[n]
slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
ax, labels = hma.plot_spectrum(n)
plot_peakslices(x_keVs, slices, ax=ax);
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')
peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]
multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel);
Mm, these histogram equalized peak maps aren't very informative. Most of the slices are just noise. I can not be sure how manganese and iron are correlated. Need to plot that in a different way.
FeKa_map = peak_maps[3]
MnKa_map = peak_maps[0]
fig, [ax, ax1] = plt.subplots(ncols=2, sharex=True, sharey=True, figsize=[9, 5])
ax.imshow(MnKa_map)
ax.set_title('Mn_Ka')
y, x, z = hot_pixel
ax.scatter([x], [y], c='r', marker='+')
ax1.imshow(FeKa_map, vmax=2)
ax1.set_title('Fe_Ka')
ax1.scatter([x], [y], c='r', marker='+');
My first question is, are there any other manganese speckles? And second, do the correlate with iron speckles, as I would expect?
This question needs to wait. Need speckle segmentation functions.
n = 10 # Fe_Ka
x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n]
hot_pixel =hma.hotmax_pixels[n]
slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
ax, labels = hma.plot_spectrum(n)
plot_peakslices(x_keVs, slices, ax=ax);
ax.set_ylim(-10, 1.15*y_hot.max())
plot_ptrn('Fe', -2, ax);
plot_ptrn('Ca', -1, ax);
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')
peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]
multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel);
peak_idxs = hma.peak_idxs_list[10]
x_keVs[peak_idxs[0]] - x_keVs[peak_idxs[6]]
Again we see an escape peak for Fe_Ka.
n = 3
x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n]
hot_pixel =hma.hotmax_pixels[n]
slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
We need to make sure we do not confuse Cl_Ka with a Rhodium L peak.
import moseley as mos
fig, [ax, ax1] = plt.subplots(nrows=2, sharex=True, figsize=[9, 4])
hma.plot_spectrum(n, ax=ax)
plot_peakslices(x_keVs, slices, ax=ax);
ax1.set_ylim([0, 1.7])
mos.XFluo('Cl', tube_keV=30).plot(ax=ax1, color='b', peak_labels='full')
mos.XFluo('Rh', tube_keV=30).plot(ax=ax1, color='r', peak_labels='full')
mos.XFluo('K', tube_keV=30).plot(ax=ax1, color='g', peak_labels='simple')
This is not the case.
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')
peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]
multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel);
Seems to be a single grain of salt?
Run into error for n=16
n = 18 # Pb
x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n]
hot_pixel =hma.hotmax_pixels[n]
slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
import moseley as mos
fig, [ax, ax1] = plt.subplots(nrows=2, sharex=True, figsize=[9, 4])
hma.plot_spectrum(n, ax=ax)
plot_peakslices(x_keVs, slices, ax=ax);
ax1.set_ylim([0, 1.7])
mos.XFluo('Pb', tube_keV=30).plot(ax=ax1, color='k', peak_labels='full')
mos.XFluo('S', tube_keV=30).plot(ax=ax1, color='r', peak_labels='full')
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')
peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]
multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel);
This is most interesting. Let's take a look at the keV map for the overlapping band of sulfur and lead. For sulfur this is the S_KL3 band at 2.311 keV. For lead this is the Pb_M5O2 band and 2.372 keV.
The best map for sulfur without lead is actually the potassium map #4
mid_keV = (2.311 + 2.371) / 2
SKa_keV_map = keV_maps[4]
fig, [ax, ax1] = plt.subplots(ncols=2, sharex=True, sharey=True, figsize=[9, 5], constrained_layout=True)
ax.imshow(peak_maps[4])
ax.set_title('Overlapping band peak_maps[4]')
pos = ax1.imshow(SKa_keV_map, vmin=2.2, vmax=2.45)
fig.colorbar(pos, ax=ax1, shrink=0.8)
ax1.set_title('Energy of maximum \nS=2.311 keV vs Pb=2.372 keV');
Better use best lead free sulfur map in potassium (n=4)?
Ok, this is the new Gaussian code base. For now we do not return the fancy peak shape parameters to deconvolve overlapping bands. Soon this will be necessary. For now let's simply take a look at the peak maps...